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ABSTRACT 

A finite difference model for elastic waves is introduced. The model is 
based on the first order system of equations for the velocities and 
stresses. The differencing is fourth order accurate on the spatial 
derivatives and second order accurate in time. The model is tested on a 
series of examples including the Lamb problem, scattering from a plane 
interfaces and scattering from a fluid-elastic interface. The scheme is shown 
to be effective for these problems. The accuracy and stability is insensitive 
to the Poisson ratio. For the class of problems considered here we find that 
the fourth order scheme requires from two-thirds to one-half the resolution of 
a typical second order scheme to give comparable accuracy. 
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Introduction 


In this paper we introduce and validate a fourth order accu- 
rate finite difference scheme for the computation of waves in an elas- 
tic medium. The method is based on the first order system for the 
linearized velocities and stresses. The numerical scheme is a fourth 
order accurate version of the MacCormack scheme (see Gottlieb and 
Turkel, 1976) which has been shown to be effective for acoustic wave 
propagation (Maestrello et. al., 1981). 

In this paper we consider only the heterogeneous formulation, 
where variable, possibly discontinuous, elastic parameters are used 
and the equation is differenced over the entire computational domain. 
For a discussion of the differences between the heterogeneous formu- 
lation and the homogeneous formulation in which interface conditions 
are explicitly imposed we refer to (Kelly et. al., 1976). In (Kelly 
et. al., 1976) a second order accurate scheme for the coupled system 
of wave equations obtained from the displacement formulation of elas- 
todynamics is introduced. We refer to this as the (2-2) scheme since 
it is second order accurate in time and second order accurate in 
space. This scheme is commonly used in seismological applications and 
we will directly compare the proposed fourth order scheme with the 
(2-2) scheme. The proposed scheme is second order accurate in time 
and fourth order accurate in space. It will be referred to as a (2-4) 
scheme. 

The first order system has been used previously to compute 
elastic waves. The reader is referred to (Clifton, 1967), (Madariaga, 
1976), (Virieux and Madriaga, 1982) and (Emerman et. al., 1982) for a 
discussion of some schemes that have been used previously. All of the 
above schemes are only second order accurate in space and are there- 
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fore distinctly different from the scheme proposed here. It is possi- 
ble to extend the (2-2) scheme of (Kelly et. al., 1976) to make the 
spatial differences fourth order accurate. This is similiar to the 
derivation of fourth order accurate schemes for the wave equation (see 
for example (Alford et. al., 1974)). Such a fourth order accurate 
scheme has been used for acoustic and S-H wave propagation (see for 
example (Frankel and Clayton, 1984)). The extension to the system of 
coupled wave equations for an elastic medium is somewhat more compli- 
cated, due to the presence of the mixed derivatives but in practice 
would be relatively straightforward. 

As far as we are aware, no such (2-4) scheme has been studied 
for elastic wave propagation. A major source of difficulty is in 
obtaining a stable implementation of the free surface condition which 
is of sufficient accuracy to maintain the overall fourth order accura- 
cy of the scheme ( see (Oliger, 1976) for a discussion of the impor- 
tance of the numerical approximation of boundary conditions). The 
need to approximate one-sided derivatives for the free surface condi- 
ton is a major numerical difference between elastodynamics and acous- 
tics where the pressure release boundary conditon does not involve 
spatial derivatives. The scheme proposed here, for which the stresses 
are dependent variables, permits a stable and accurate implementation 
of the free surface condition even for large Poisson ratios (see below 
). 

A higher order accurate scheme (in fact infinite order accu- 
rate ) has been proposed by (Kosloff et. al., 1984). This method is 
based on the Fourier pseudo-spectral method and requires a periodic 
extension of the computational domain. The comparison between this 
method and the (2-4) scheme proposed here remains to be carried out. 
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in particular for realistic size seismological models involving dis- 
continuous elastic parameters. 

The proposed (2-4) scheme is based on operator splitting whe- 
re the two-dimensional equations are solved first in one dimension and 
then in the second. The concept of operator splitting and the details 
of the (2-4) scheme are discussed in the next section. We then pres- 
ent numerical examples for the Lamb problem, scattering from a hori- 
zontal interface and scattering from a fluid-elastic interface. We 
finally present our conclusions. 
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NUMERICAL SCHEME 

This section is divided into four parts. In part A we dis- 
cuss the formulation of linear elastodynamics as a first order system 
for the velocities and stresses. In part B, the numerical scheme is 
described. In part C, we consider the choice of timestep and in part 
D the boundary conditions are described. 


A. Formulation 

The equations of linear isotropic elastodynamics in Cartesian 
coordinates are 


pu t = 

r ll,x + 

T i2,y 

pv t = 

r 12 , x + 

T 22 , y 

r ll,t 

= ( A + 

2p )u x 

r 22 , t 

= pu x 

+ ( x 

T 12 , t 

= p(v x 

* V 


( 1 ) 


In equation (1) u and v are the horizontal and vertical velocities 
respectively, are the components of stress tensor. The elastic 

parameters are the density p and the Lame constants A and u. The com- 
pressional and shear wave speeds Cp and C g are given by 


(A + 2 u) 
P 


( 2 ) 


_ £ 
P 


and C > C , 
P s 


Associated with these speeds there are the spatial wave- 
lengths ^p = C p /f, A s =C g /f where f is a characteristic frequency. Hence 
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JB = 


k >u 


The Poisson ratio is given by 


v 


1 

2 


1 - 2C^/C 

1 “ C s /C n 
s p 


L 1 


so that 



2(l-») 
1-2 * ' 


(3) 


(4) 


(5) 


For seismic problems it is frequently assumed that v is 

around 0.25 and thus that C /C =/3~. It then follows that the shear 

P s 

wavelengths are about 60% smaller than the compress ional wavelengths. 
Therefore the spatial resolution requirements must be based on the 
shorter shear wavelengths. In addition there exist interface and sur- 
face waves which generally have smaller wavelengths than the shear 
waves. There are applications where these waves are not generally 
considered to be important. Nevertheless, it is necessary to resolve 
these waves sufficiently well to prevent a spurious transfer of energy 
into other wave modes. We also do not want numerical dispersion for 
these waves to interfere with the generation and interpretation of the 
waves of interest. In some cases, e.g. weathering layers, the Poisson 
ratio can be considerably larger that 0.25. This further accentuates 
the difference in resolution requirements between shear and compres- 
sional waves . 

For convenience we write equation (1) in vector form 

W. = ’ AW + BW t , (6) 

z x y 

T • 

11' t 22 ' T 12' and the matr ^ ces A and B are 9 lven by 


where W = (u, v, t 
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We can also write equation (6) as a symmetric system 


where 


E 0 W t 


A 0 W x + B 0 W y 




0 

P 

0 

0 

0 


0 

0 

X+2 u 

(X+2y)^— 

(\+2u) 2 -u 2 

0 


0 

0 

~u 

( \+2^ ) 2 —u 2 ’ 
X+2 ii 

( \+2u ) 2 — y 2 
0 




0 

0 

0 

0 

1 


1 

0 

0 

0 

0 


0 

0 

0 

0 

0 


(7) 
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We note that the matrices E qAq and Eq^Bq do not commute. 
Since p, X, p are independent of time, Eq is also independent of time. 
The matrices Aq and Bq are also independent of the spatial 
coordinates. It follows that equation (7) can be be written in the 
general form 

u t ■ £ x * ”y < 81 

where the vector U is EqW and the functions f and h are AqW and BqW 
respectively. Equations of the form of (8) are called 
divergence-free. We will develop a numerical scheme for equations of 
this form. Since equation (7) was obtained from equation (6) by line- 
ar manipulations the resulting scheme will be valid for equation (6). 
The numerical code is based on equation (6). 

B. Numerical Algorithm 

The numerical algorithm will be based on the method of dimen- 
sional splitting (Strang, 1968). In this method, the two-dimensional 
equation (6) is updated for one timestep by first solving the equation 
in the x-direction and then in the y direction. For the next timestep 
the order of the x and y updates is reversed. In order to advance the 
solution from time level n to level n+2 we use the formula 

w " t2 ■ L * L y L y L * wn; (9) 


L x ,Ly represent the solution to the one dimensional problems 
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W. = AW 
t x 

W. = BW < 1Q ) 

t y 

respectively. 

It follows from equation (8) that to advance the solution one 
timestep we solve equation (6) using only the terms involving the 
x-der ivatives . We then use these new values for the dependent vari- 
ables and solve equation (6) using only the terms involving the 
y-derivatives. This gives the solution at time level n+1. To update 
to the next time level we repeat the procedure but reverse the order 
of the x and y updates. By reversing the order the resulting scheme 
is second order accurate in time. 

The use of operator splitting has the following advantages. 

(a) The maximum allowed timestep is governed by 

the associated one-dimensional equations. 

The allowed timestep is larger than 
for unsplit schemes. 

(b) split schemes are particularly well suited to 

multi-processing (see Schenck et. al.,1985)). 
For example in the first step all 
the x updates for each line y = constant can 
be done independently. 

(c) split schemes tend to have smaller phase 

errors than schemes (Turkel, 1974). 

The proposed scheme can be implemented without splitting if desired 
(see (Gottlieb and Turkel, 1976)). 

We need to specify a one-dimensional scheme for the operators 
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L and L . We use a scheme that is second order in time and fourth 
x y 

order in space (Gottlieb and Turkel, 1976). For the equation 


n = 

f 

(11) 

t 

X 


Cl 

*-*• 

II 

(7(£ ? + i - - “U - f w>> 

(12) 

< 

■ 2 (U i +U i + Mx (7(f i f i-l> - (f i-l 

- fi_ 2 )>) 


alternating with a symmetric variant. In equation (12) the subscript 
' i * denotes the spatial grid point and the superscript 'n' denotes the 
time level. The resulting scheme is stable when applied to equation 
(6) provided 

C < .67 . (13) 

Ax p 

The ratio on the left hand side of equation (13) will be called the 
CFL number. 

The scheme described by equation (12) together with operator 
splitting has been applied successfully to a variety of wave propa- 
gation problems. The reader is referred to (Maestrello et. al. r 1981) 
for an application to the computation of acoustic disturbances in a 
jet. As far as we are aware, split schemes have not yet been applied 
to elasto-dynamic problems. The scheme based on equation (12) is dis- 
sipative and damps the high frequencies (Gottlieb and Turkel, 1976). 

It has a greatly reduced truncation error compared to second order 
schemes . 

The order of accuracy in time is second order. Hence, one 

2 

achieves true fourth order accuracy only if At = 0(Ax ). Thus in gen- 
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eral the timestep must be smaller than the maximum timestep allowed by 
stability (equation (13)) in order to prevent the errors due to the 
second order accuracy in time from dominating the truncation error. 

In elastic wave propagation, however, it is possible to choose a 
timestep near that allowed by stability. This is due to the differ- 
ence in size between the compressional and shear velocities. We have 
verified this by extensive numerical experiments and provide a justi- 
fication of this immediately below. 


C. Choice of Timestep 

In elasticity the spatial resolution and the timestep are 
chosen according to two different considerations. This is because of 
the existence of two different characteristic speeds. In general, the 
spatial resolution is based on the shear velocity while the stability 
restriction , equation (13), is based on the larger compressional 
velocity. This allows larger timesteps without losing temporal accu- 
racy. 

In order to see this let f be a characteristic frequency of 
the problem so that the associated compressional and shear wavelengths 
satisfy equation (3). Assume that we could completely decouple the 
compressional and shear waves and use a different mesh size for each 
wave type that is appropriate to obtain a specified accuracy. The 
timesteps for each wave type would then be determined by 


At 

AX P 

At 



(14) 


where K is the CFL number. For the scheme described in equation (12) 
the maximum value of K allowed by stability is 0.67. In order to 
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reduce temporal errors by reducing the timestep, K is reduced. It 
follows from equation (14) that 


Ax 

Ax 


3b 

C s 


(15) 


This is simply a statement that a smaller mesh is required to resolve 
shear waves. In practice we must use the same mesh, Ax g , for both 
compress ional and shear waves. Since the timestep for the system is 
based only on Cp according to equation (13) we have 

K' Ax 

At = ? a (16) 

P 

where K' is the CFL number used in the actual computation. Comparing 
equation (14) with equation (16) we find that 

KC _ 

K' = p- 2 . (17) 

s 

According to equation (17) the CFL number that is actually used in the 
code (K' ) is equal to the CFL number based on temporal accuracy (K) 
multiplied by C /C (~/3). Hence we should choose the timestep to be 
about 60% larger than the timestep that would give acceptable results 
for the compressional and shear waves decoupled. It is clear from 
this argument that the compressional waves will generally be oversam- 
pled but this is unavoidable for any scheme that solves the equation 
( 1 ). 

It follows from the discussion above that the spatial mesh 
size is governed by the shear velocity while the explicit timestep 
restriction due to stability is governed by the larger compressional 
velocity. This is essentially the case for all explicit finite dif- 
ference schemes. For a scheme that is second order accurate in both 
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space and time, the optimal timestep is usually the maximum allowed by 
stability. For a scheme that is second order accurate in time and 
fourth order accurate in space the optimum time step is typically 
smaller than that allowed by stability. This is in order to prevent 
the second order temporal errors from dominating the truncation error. 
However when the explicit timestep restriction itself forces a small 
timestep, the (2-4) scheme can be run at or near its stability limit 
and the temporal errors will be negligible. The use of operator 
splitting then permits a larger timestep as indicated above. 

As the Poisson ratio increases the ratio of the compressional 
velocity to the shear velocity becomes larger. The corresponding 
timestep forced by the stability restriction can become much smaller 
than that required for accuracy. We have found that for a Poisson 
ratio around 0.3 the (2-4) scheme can be run at its stability limit 
(equation (13)) with little degradation in accuracy. For very high 
Poisson ratios, implicit schemes such as that proposed by (Emerman et. 
al., 1982) would become more efficient since they would permit times- 
teps that are not restricted by stability. In practice this must be 
balanced against the additional work per timestep required by implicit 
schemes. This discussion highlights a difference between elasticity 
and acoustics or S-H wave propagation that makes (2-4) schemes more 
efficient in the elastic case. The discussion would be applicable to 
acoustic wave propagation only if the computational model included 
wide variations in the size of the acoustic velocity. 

The numerical experiments using this (2-4) scheme were per- 
formed in a rectangular region 

OSxSI^; -L 2 ^yS0. (18) 

A Cartesian mesh with Ax=Ay is always used independent of the shape of 
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the interfaces. It requires 112 floating point operations per grid 
point (multiplications and adds) to advance the solution from t to 
t+At (ignoring boundary calculations). The simplicity of the algo- 
rithm and the mesh makes it easy to treat many different applications. 
Changes in geometry or even in the equations do not introduce any sig- 
nificant complications. 


D. Boundary Conditions 

The interior scheme, given by equation (12) requires the use 
of two points in every direction from the point being advanced. At 
boundaries some numerical boundary treatment must be introduced 
because the boundary point does not have enough neighbors to implement 
the scheme. In addition, there are physical boundary conditions that 
must be imposed. 

In each sweep all boundary fluxes are defined at points out- 
side of the computational domain by a third order extrapolation. For 
example, for the scheme defined by equation (12) if i=n is a boundary 
point we define, 


(19) 


The extrapolated fluxes f n+1 and f n+2 are used to complete the forward 
predictor. Other boundaries are treated similarly. The third order 
extrapolation given by equation (17) enables us to compute a solution 
at the boundary which maintains the overall fourth order accuracy of 
the scheme. It is shown in (Oliger, 1976) that the accuracy obtained 
from fourth order accurate interior schemes can be significantly 
degraded if a third order accurate boundary treatment is not used. 


n+1 

= 4f 

— 6f , 

+ 4f 0 

- f ~ 

n 

n-1 

n-2 

n— 3 

n+2 

= 4f n + l " 6f n 

+ 4f n-l 

- f n— 2 
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It remains to impose the correct physical boundary conditions 
which are necessary to complete the specification of the ini- 
tial boundary value problem given by equation (1) in the region given 
by equation (18). In typical calculations there are two types of 
boundaries. The upper boundary, y=0, is a free surface along which 
the appropriate conditions are T i? =T 22 = ® ' Since the scheme is a cen- 
tered scheme we must calculate the other variables u, v, by some 

other means. We calculate these variables by extrapolation based on 
characteristic variables. 

Based on a one dimensional analysis (i.e. neglecting x deriv- 
atives) the quantity 

= t/p/i u + r ^2 (20) 

is convected toward the boundary with velocity C s , 

R 2 = v'(X + 2 p)p v + T 22 (21) 

is convected toward the boundary with velocity Cp and 

R 3 = (X + 2v)t 1 ^-ht 22 ( 22 ) 

moves with zero normal velocity. We specify R^ , R 2 , Rg by the values 
obtained from the interior scheme with the boundary fluxes extrapo- 
lated using equation (19). This combined with the free surface condi- 
tions t 22 =t 12 =0 enables us to obtain all of the dependent variables. 
(For problems with a surface source, the above procedure is carried 
out with t 22 and T i 2 ec I ual to the prescribed forces). The use of 
one-dimensional characteristic variables was shown in (Gottlieb et. 
al., 1982) to enhance the stability of the boundary treatment, we 
have found that the specification of the zero characteristic R 3 is 
necessary when the Poisson ratio is large or when the material is 
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highly anisotropic. 

In addition to the free surface there are three other bounda- 
ries; x=0, x=L 1 and y=-L 2< These are all artificial boundaries and 
appropriate absorbing boundary condition must be imposed. We use the 
one-dimensional characteristics given by equations (20-22) (or the 
analogue obtained by neglecting y derivatives when the boundary is the 
line x = constant ) and impose the boundary condition that the incom- 
ing characteristic variable is zero. This characteristic condition is 
the same as the viscous boundary condition proposed by (Lysmer and 
Kuhlemeyer, 1969). Using the first order system (equation (1)) it 
becomes a Dirichlet condition which can be implemented very easily. 
These conditions are exact for compressional and shear waves impirrging 
normally on the boundary. Since they are based on one-dimensional 
characteristics they do not absorb Rayleigh waves. The Rayleigh wave 
reflection can be ameliorated by burying the source or by making a 
subsidiary calculation of the free space problem with just the upper 
layer of the model and subtracting that solution from the computed 
solution. We are investigating more satisfactory solutions to this 
problem. 


In (Clayton and Engquist, 1977) boundary conditions based on 
paraxial approximations to elastic waves were introduced. The first 
order condition in (Clayton and Engquist, 1977) is equivalent to the 
characteristic condition in one dimension but differs in two dimension 
because of the replacement of the shear stress by a displacement. The 
higher order paraxial boundary conditions are theoretically more 
effective absorbers as the waves deviate from normal incidence. (Co- 
hen et. al., 1981) have found that these higher order boundary condi- 
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tions are not significantly more effective than the viscous boundary 
conditions and we have not found pernicious reflections in the exam- 
ples considered in the next section. None of these boundary condi- 
tions accounts for the cylindrical decay of the compressional and 
shear waves. A family of boundary conditions which account for the 
cylindrical decay of acoustic waves is described in (Bayliss and Turk- 
el, 1980). 
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NUMERICAL RESULTS 

In this section we present some numerical examples which 
illustrate the advantages of the finite difference scheme described 
above. In the first example we consider the Lamb problem for a sur- 
face line source. Specifically we take 

r 22 = f ( t ) 6 (x) (23) 

on the free surface. The functional form of the source is 

. —2 ( t— t ) 2 / a 2 — 2 ( t /a ) ) 2 

f(t) = - (e s -e S ) (t-t s ).(24) 

a 

The material properties were chosen so that the compressional speed 
was 7000 ft. /sec., the shear speed was 4000 ft. /sec. and the density 
was taken as 1.0. The parameters of the source were a = .017 sec and 
the time shift (t g ) was 0.285 sec. The spatial delta function was 
modelled by the discrete function which is zero except at the location 
of the source where it is equal to 1/h (h is the grid size). 

The coupled wave equation for the elastic displacements can 
also be considered as coupled wave equations for the elastic veloci- 
ties provided we use a source that is the derivative of equation (24). 
Using the parameters described above the peak frequency of the deriva- 
tive of the source is approximately 26.5 Hertz. We compared the sol- 
ution with the exact solution (Miklowitz, 1978) and with solutions 
obtained from the (2-2) scheme (Kelly et. al., 1976). A second order 
accurate implementation of the free surface condtion was used (B. Nair 
private communication). 

In Figures la and lb we compare the horizontal velocity u at 
a receiver location of 100 ft. from the source. We compare the exact 
solution with results from the finite difference calculations using a 
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mesh size of 10 ft. It is apparent that the (2-4) scheme is more 
accurate. In refining the grid we found that the (2-2) scheme with 5 
ft. spacing gave a solution that was comparable to the (2-4) scheme 
with 10 ft. spacing. This is exhibited in Figure lc. 

In Figure 2 we demonstrate the effect of increasing the Pois- 
son ratio. The parameters are the same except that the compressional 
velocity is increased to 16000 ft. /sec. giving a Poisson ratio of 
0.4667. The results demonstrate that the (2-4) scheme is stable and 
that there is little loss in accuracy as the Poisson ratio increases. 
The implementation of the free surface condition which we used with 
the (2-2) scheme was not stable at this Poisson ratio. We did not 
test the implementation of (I lan and Loewenthal, 1976) which appears 
to be more stable at higher Poisson ratios. 

We next demonstrate the accuracy of the numerical scheme on 
body waves reflected from an interface. We consider a surface source 
of the above form and an interface 400 ft. below the free surface. At 
the interface the velocities jump to 12000 ft. /sec. (compressional) 
and 8000 ft. /sec. (shear). The density is kept constant. 

We have found that for a given source and mesh size the Ray- 
leigh waves tend to be significantly less accurate than the body 
waves. There are several reasons for this. The Rayleigh waves tend 
to be a higher order derivative than the body waves (see Miklowitz, 
1978) and thus there is more energy in the higher frequencies. In 
addition, the Rayleigh wave profiles are steep in the normal direction 
thus requiring more resolution. The accuracy is also sensitive to 
the treatment of the free surface condition. Finally Rayleigh waves 
propagate parallel to the grid and many schemes tend to be least accu- 
rate for waves that propagate parallel to the grid. 
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For some problems in seismology the reflected body waves are 

more important than the Rayleigh waves. In order to asses 
•the accuracy of the two schemes on the body waves we compared sol- 
utions on the free surface obtained with both numerical schemes. The 
time intervals were chosen so as not to include the Rayleigh wave. 

Mode converted arrivals and multiple reflections off the free surface 
did occur during the chosen time interval. In order to enhance the 
separation of the arrivals, the duration of the source pulse was 
reduced by setting a = .0125 sec. 

The programs were run with grid refinement until no further 
changes were visible. It was found that the (2-4) scheme with a mesh 
size of 5 ft. could be taken as the exact solution. A comparision for 
the horizontal velocity u at a receiver location of 100 ft. is shown 
in Figures 3a and 3b. In Figure 3a we compare the solution obtained 
with the (2-4) scheme with a mesh size of 5 ft. to the solution 
obtained with the (2-4) scheme and a mesh size of 7.5 ft. It is 
apparent that the solutions agree closely, with the major source of 
errors being a slight amplitude attenuation in the later arrivals. In 
contrast, we compare in Figure 3b the (2-2) scheme with a mesh size of 
7.5 ft. to the (2-4) scheme with a mesh size of 5.0 ft. The second 
order solution is quite inaccurate, exhibiting both amplitude errors 
and additional lobes characteristic of numerical dispersion. 

In order to assess the properties of the two schemes for 
waves propagating at wider angles to the grid we make the same compar- 
ison in Figures 4a and 4b for a receiver location at 800 ft. The same 
conclusion can be drawn. In Figure 4c we compare the (2-4) scheme and 
the (2-2) scheme both with mesh sizes of 5 ft. It can be inferred 
that the second order scheme with a mesh size of 5 ft. is comparable 
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in accuracy to the fourth order scheme with a mesh size of 7.5 ft. 

In summary the results for this model demonstrate that both 
schemes perform worst for waves propagating nearly parallel to the 
grid. The (2-4) scheme is less anisotropic. For models of this size 
and sources of the form of equation (24) body waves require roughly 
two-thirds less resolution in each dimension with the (2-4) scheme 
than with the (2-2) scheme for comparable accuracy. Rayleigh waves 
require roughly half the resolution. The accuracy and stability of 
the (2-4) scheme does not appear to be sensitive to Poisson ratio. We 
have obtained results for problems with curved interfaces, including 
models which allow caustic formation. The conclusions for these cases 
are similiar to those already presented. 

In our final example we consider the problem of Rayleigh wave 
scattering from a fluid. The geometry is shown in Figure 5. For sim- 
plicity, the Rayleigh wave is generated by a surface source in the 
elastic region. The acoustic fluid is treated by setting the shear 
modulus to zero and differencing across the interface. The free sur- 
face condition is replaced by a pressure release boundary condition 
over the fluid. 

This problem is designed to demonstrate that the scheme is 
accurate and stable at a fluid-elastic interface. We refer to (Steph- 
en, 1983) for a discussion of the numerical difficulties that can be 
expected at a fluid-elastic interface. The source function is the 
same as for the first example. In Figure 6 we plot a snapshot of the 
vertical velocity at t = 0.35 sec. The impinging Rayleigh wave gener- 
ates an interfacial wave which travels along the fluid-elastic inter- 
face and a wave which travels throught the fluid. The slow velocity 
of the interfacial wave governs the resolution requirements of the 
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problem. Body waves are generated at the interface and propagate away 
from the fluid. 

In order to assess the accuracy of the computation, we exam- 
ine time traces of the horizontal velocity for different grids. We 
choose a receiver location at 1500 ft. so that interference by the 
slow interfacial wave can be seen. In Figures 7a and 7b the solution 
with grid sizes of 4 and 6 ft. are compared with the solution obtained 
with a grid size of 2.5 ft. The solution does not change if the grid 
is refined further. It is apparent that the solution converges as the 
grid is refined. We have not run this case with the (2-2) scheme but 
based on the results presented above we expect that the (2-2) scheme 
would require substantially more resolution to obtain a solution of 
comparable accuracy. Considerably more resolution is required than 
for the corresponding Lamb problem. This is primarily because of the 
slow velocity of the interfacial wave. The computations were stable 
provided the CFL number (see above ) was reduced to 0.5. 

This computation, together with others that we have made, 
demonstrates that there are no stability problems in applying the 
scheme at a fluid-elastic interface. The change to a pressure release 
boundary condition over the fluid provides a further test of the 
robustness of the scheme. We have repeated this computation with the 
Poisson ratio in the elastic region increased to 0.4667 and again no 
instabilities were found. Although the resolution requirements in 
this problem are governed by a slow interfacial wave, we have examined 
body waves reflected from a fluid-elastic interface and found no 
degradation in accuracy compared to a corresponding elastic-elastic 
interface. The numerical scheme provides a simple and robust method 
to compute scattering from curved fluid-elastic interfaces. 
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CONCLUSION 

The numerical results presented here clearly demonstrate the 
improvements in accuracy that can be obtained from the use of fourth 
order accurate differences to compute elastic waves. We have tested 
the scheme on a variety of problems involving wave scattering from 
curved interfaces, fluid-elastic interfaces and regions of high Pois- 
son ratio. The results presented here are representative of the 
improvements due to the fourth order differencing. 

The scheme presented here permits an accurate and robust 
implementation of the free surface condition. For body waves, with 
sources of the same form as equation (24), we have found that from 
13-16 points per shear wavelength, based on the peak frequency of the 
derivative of the source, is required for reasonable accuracy. This 
is for models approximately 20-30 wavelengths in size. Larger models 
require more resolution per wavelength. This is shown mathematically 
by (Kreiss and Oliger, 1973) and demonstrated numerically by (Stephen, 
1982). Resolution requirements grow at a much slower rate with fourth 
order differencing (Kreiss and Oliger, 1973). Considerably more 
resolution is required to approximate Rayleigh waves than body waves. 

The scheme is suitable for vector computers. A sustained 
rate of 56 MFLOPS (Million Floating Point Operations Per Second) can 
be obtained on a Cray 1-S. On a Cray XMP/48 using only one processor 
a rate of 83 MFLOPS can be obtained. Using all four processors rather 
than one a sustained rate of 328 MFLOPS can be obtained. The improve- 
ment in using four processors rather than one is 3.95. (The authors 
are grateful to Michael Booth of Cray Research Inc. for assistance in 
adapting the program to run on multi-processors). Schemes based on 
operator splitting are particularly well suited to multi-processing. 
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For example during an x update each line y=constant can be given to a 
separate processor. A further discussion of this can be found in 
(Schenck et. al., 1985). A benchmark problem with a grid of 701x301 
and run for 2560 timesteps, required 660 seconds on the Cray 1-S and 
116 seconds on the XMP/48. (This run and all of the other cases dis- 
cussed here used an expanding grid as suggested by (Boore, 1972)). 

There are applications where the explicit imposition of 
interface conditions may be more efficient than the heterogeneous for 
mulation. For example, models with localized regions of low velocity 
would require significant oversampling unless different grids were 
used in different regions. The scheme described here appears to be 
applicable to a wide variety of problems of seismological interest. 

It could also be used with the explicit imposition of interface condi 
tions (homogeneous formulation) provided interface conditions of suf- 
ficient accuracy were imposed. 
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FIGURE CAPTIONS 


l.a Comparison of exact solution to the Lamb problem (solid curve) 
to solution obtained with (2-4) scheme using a grid size of 
10 ft. (dotted curve). Horizontal velocity (u) plotted at a 
receiver location 100 ft. from the source. 

l.b Comparison of exact solution to the Lamb problem (solid curve) 
to solution obtained with (2-2) scheme using a grid size of 
10 ft. (dotted curve). Horizontal velocity (u) plotted at a 
receiver location 100 ft. from the source. 

l.c Comparison of exact solution to the Lamb problem (solid curve) 
to solution obtained with (2-2) scheme using a grid size of 
5 ft. (dotted curve). Horizontal velocity (u) plotted at a 
receiver location 100 ft. from the source. 

2 Comparison of exact solution to the Lamb problem (solid curve) 
to solution obtained with (2-4) scheme using a grid size of 
10 ft. (dotted curve). Horizontal velocity (u) plotted at a 
receiver location 100 ft. from the source. Poisson ratio is 
0.4667. 


3. a Comparison of reflected body waves. Solid curve is (2-4) scheme 
using a mesh size of 5 ft. Dotted curve is obtained from (2-4) 
scheme using a mesh size of 7.5 ft. Horizontal velocity (u) 
plotted at a receiver location 100 ft. from the source. 
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3. b Comparison of reflected body waves. Solid curve is (2-4) scheme 

using a mesh size of 5 ft. Dotted curve is obtained from (2-2) 
scheme using a mesh size of 7.5 ft. Horizontal velocity (u) 
plotted at a receiver location 100 ft. from the source. 

4. a Comparison of reflected body waves. Solid curve is (2-4) scheme 

using a mesh size of 5 ft. Dotted curve is obtained from (2-4) 
scheme using a mesh size of 7.5 ft. Horizontal velocity (u) 
plotted at a receiver location 800 ft. from the source. 

4.b Comparison of reflected body waves. Solid curve is (2-4) scheme 
using a mesh size of 5 ft. Dotted curve is obtained from (2-2) 
scheme using a mesh size of 7.5 ft. Horizontal velocity (u) 
plotted at a receiver location 800 ft. from the source. 

4.c Comparison of reflected body waves. Solid curve is (2-4) scheme 
using a mesh size of 5 ft. Dotted curve is obtained from (2-2) 
scheme using a mesh size of 5 ft. Horizontal velocity (u) 
plotted at a receiver location 800 ft. from the source. 

5 Model for problem of Rayleigh wave scattering from fluid. 

6 Snapshot of vertical velocity (v) obtained for problem of 
Rayleigh wave scattered from fluid. Time is 0.35 sec 

7. a Comparison of solutions obtained for problem of Rayleigh wave 
scattering from fluid. Solid curve is obtained with grid size 
of 2.5 ft. Dotted curve had grid size of 6 ft. Horizontal 



29 


velocity (u) plotted at receiver location 1500 ft. from source 

7.b Comparison of solutions obtained for problem of Rayleigh wave 
scattering from fluid. Solid curve is obtained with grid size 
of 2.5 ft. Dotted curve had grid size of 4 ft. Horizontal 
velocity (u) plotted at receiver location 1500 ft. from source 
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